function banded_systems()
% BANDED_SYSTEMS 带状系统求解演示
% 
% 演示三对角、五对角等带状矩阵的高效求解方法

fprintf('=== 带状系统求解 ===\n\n');

%% 1. 三对角系统
fprintf('1. 三对角系统\n');
fprintf('三对角矩阵只有主对角线和相邻的两条对角线非零\n\n');

% 创建三对角矩阵
n = 5;
a = -ones(n-1, 1);  % 下对角线
b = 4*ones(n, 1);   % 主对角线
c = -ones(n-1, 1);  % 上对角线

A_tri = diag(a, -1) + diag(b, 0) + diag(c, 1);
fprintf('三对角矩阵 A =\n');
disp(A_tri);

% 右端向量
b_rhs = [1; 2; 3; 4; 5];
fprintf('右端向量 b = [%g; %g; %g; %g; %g]\n', b_rhs);

% Thomas算法求解
x_thomas = thomas_algorithm(a, b, c, b_rhs);
fprintf('Thomas算法解: x = [%.6f; %.6f; %.6f; %.6f; %.6f]\n', x_thomas);

% 验证解
residual = A_tri * x_thomas - b_rhs;
fprintf('残差 ||A*x - b|| = %e\n', norm(residual));

% 与MATLAB内置求解器比较
x_matlab = A_tri \ b_rhs;
fprintf('MATLAB解差异: %e\n', norm(x_thomas - x_matlab));

%% 2. Thomas算法的稳定性
fprintf('\n2. Thomas算法的稳定性\n');
fprintf('Thomas算法在对角占优时是稳定的\n\n');

% 对角占优的三对角矩阵
fprintf('对角占优情况:\n');
a_stable = -ones(n-1, 1);
b_stable = 3*ones(n, 1);  % |b_i| > |a_i| + |c_i|
c_stable = -ones(n-1, 1);

A_stable = diag(a_stable, -1) + diag(b_stable, 0) + diag(c_stable, 1);
fprintf('对角占优矩阵:\n');
disp(A_stable);

% 检查对角占优条件
fprintf('对角占优检查:\n');
for i = 1:n
    row_sum = 0;
    if i > 1, row_sum = row_sum + abs(a_stable(i-1)); end
    if i < n, row_sum = row_sum + abs(c_stable(i)); end
    
    fprintf('行%d: |%g| > %g ? %s\n', i, b_stable(i), row_sum, ...
            mat2str(abs(b_stable(i)) > row_sum));
end

% 非对角占优的情况
fprintf('\n非对角占优情况:\n');
a_unstable = -2*ones(n-1, 1);
b_unstable = ones(n, 1);  % |b_i| < |a_i| + |c_i|
c_unstable = -2*ones(n-1, 1);

fprintf('可能不稳定的系统 (主对角元较小)\n');
x_unstable = thomas_algorithm(a_unstable, b_unstable, c_unstable, b_rhs);
fprintf('解: x = [%.6f; %.6f; %.6f; %.6f; %.6f]\n', x_unstable);

%% 3. 五对角系统
fprintf('\n3. 五对角系统\n');
fprintf('五对角矩阵有5条非零对角线\n\n');

% 创建五对角矩阵
n = 6;
d1 = ones(n-2, 1);    % 第二下对角线
d2 = -2*ones(n-1, 1); % 第一下对角线
d3 = 6*ones(n, 1);    % 主对角线
d4 = -2*ones(n-1, 1); % 第一上对角线
d5 = ones(n-2, 1);    % 第二上对角线

A_penta = diag(d1, -2) + diag(d2, -1) + diag(d3, 0) + diag(d4, 1) + diag(d5, 2);
fprintf('五对角矩阵 A =\n');
disp(A_penta);

b_penta = ones(n, 1);
fprintf('右端向量 b = ones(%d, 1)\n', n);

% 使用MATLAB求解（作为参考）
x_penta_matlab = A_penta \ b_penta;
fprintf('MATLAB解: x = [%.6f; %.6f; %.6f; %.6f; %.6f; %.6f]\n', x_penta_matlab);

%% 4. 带状LU分解
fprintf('\n4. 带状LU分解\n');
fprintf('带状矩阵的LU分解保持带状结构\n\n');

% 对三对角矩阵进行LU分解
A_band = diag(-ones(n-1, 1), -1) + diag(4*ones(n, 1), 0) + diag(-ones(n-1, 1), 1);
[L_band, U_band] = lu_banded(A_band, 1, 1);  % 下带宽=1, 上带宽=1

fprintf('原三对角矩阵:\n');
disp(A_band);
fprintf('L矩阵:\n');
disp(L_band);
fprintf('U矩阵:\n');
disp(U_band);

% 验证分解
LU_product = L_band * U_band;
fprintf('LU分解误差: %e\n', norm(LU_product - A_band, 'fro'));

%% 5. 循环三对角系统
fprintf('\n5. 循环三对角系统\n');
fprintf('循环三对角矩阵在角落有额外的非零元素\n\n');

% 创建循环三对角矩阵
A_cyclic = A_tri;
A_cyclic(1, n) = -1;  % 右上角
A_cyclic(n, 1) = -1;  % 左下角

fprintf('循环三对角矩阵:\n');
disp(A_cyclic);

% 使用Sherman-Morrison公式求解
x_cyclic = solve_cyclic_tridiagonal(a, b, c, -1, -1, b_rhs);
fprintf('循环三对角解: x = [%.6f; %.6f; %.6f; %.6f; %.6f]\n', x_cyclic);

% 验证
residual_cyclic = A_cyclic * x_cyclic - b_rhs;
fprintf('残差: %e\n', norm(residual_cyclic));

%% 6. 性能比较
fprintf('\n6. 性能比较\n');
fprintf('比较不同求解方法的性能\n\n');

sizes = [100, 500, 1000, 2000];
fprintf('矩阵规模    Thomas算法    LU分解    直接求解\n');
fprintf('--------    ----------    ------    --------\n');

for i = 1:length(sizes)
    n_test = sizes(i);
    
    % 生成大规模三对角矩阵
    a_test = -ones(n_test-1, 1);
    b_test = 4*ones(n_test, 1);
    c_test = -ones(n_test-1, 1);
    A_test = diag(a_test, -1) + diag(b_test, 0) + diag(c_test, 1);
    b_test_rhs = randn(n_test, 1);
    
    % Thomas算法
    tic;
    x_thomas_test = thomas_algorithm(a_test, b_test, c_test, b_test_rhs);
    time_thomas = toc;
    
    % LU分解
    tic;
    [L_test, U_test] = lu(A_test);
    x_lu_test = U_test \ (L_test \ b_test_rhs);
    time_lu = toc;
    
    % 直接求解
    tic;
    x_direct_test = A_test \ b_test_rhs;
    time_direct = toc;
    
    fprintf('%6d      %8.4f s    %6.4f s    %6.4f s\n', ...
            n_test, time_thomas, time_lu, time_direct);
end

%% 7. 应用示例：一维热传导方程
fprintf('\n7. 应用示例：一维热传导方程\n');
fprintf('离散化一维热传导方程得到三对角系统\n\n');

% 问题设置
L = 1;          % 区域长度
n_grid = 21;    % 网格点数
h = L / (n_grid - 1);  % 网格间距
alpha = 1;      % 热扩散系数

fprintf('网格间距 h = %.4f\n', h);
fprintf('网格点数 = %d\n', n_grid);

% 构造系数矩阵（内部点）
n_inner = n_grid - 2;  % 内部点数
a_heat = -alpha * ones(n_inner-1, 1) / h^2;
b_heat = (2*alpha / h^2) * ones(n_inner, 1);
c_heat = -alpha * ones(n_inner-1, 1) / h^2;

% 右端项（热源）
f = @(x) sin(pi * x);  % 热源函数
x_grid = linspace(h, L-h, n_inner);  % 内部网格点
b_heat_rhs = f(x_grid)';

% 边界条件 u(0) = 0, u(L) = 0 已经包含在系统中

% 求解
u_inner = thomas_algorithm(a_heat, b_heat, c_heat, b_heat_rhs);

% 添加边界值
u_full = [0; u_inner; 0];
x_full = linspace(0, L, n_grid);

fprintf('数值解在中点 x=0.5 处的值: %.6f\n', u_full(n_grid/2 + 1));

% 绘制解
figure('Name', '一维热传导方程数值解');
plot(x_full, u_full, 'bo-', 'LineWidth', 2, 'MarkerSize', 6);
xlabel('位置 x');
ylabel('温度 u(x)');
title('一维热传导方程的数值解');
grid on;

% 解析解比较（如果存在）
u_exact = sin(pi * x_full) / (pi^2 * alpha);
hold on;
plot(x_full, u_exact, 'r--', 'LineWidth', 2);
legend('数值解', '解析解', 'Location', 'best');

error_max = max(abs(u_full - u_exact'));
fprintf('最大误差: %e\n', error_max);

end

function x = thomas_algorithm(a, b, c, d)
% THOMAS_ALGORITHM Thomas算法求解三对角系统
% 输入:
%   a - 下对角线 (长度 n-1)
%   b - 主对角线 (长度 n)
%   c - 上对角线 (长度 n-1)
%   d - 右端向量 (长度 n)
% 输出:
%   x - 解向量

n = length(b);
x = zeros(n, 1);

% 前向消元
c_prime = zeros(n-1, 1);
d_prime = zeros(n, 1);

c_prime(1) = c(1) / b(1);
d_prime(1) = d(1) / b(1);

for i = 2:n-1
    denominator = b(i) - a(i-1) * c_prime(i-1);
    c_prime(i) = c(i) / denominator;
    d_prime(i) = (d(i) - a(i-1) * d_prime(i-1)) / denominator;
end

d_prime(n) = (d(n) - a(n-1) * d_prime(n-1)) / (b(n) - a(n-1) * c_prime(n-1));

% 后向替换
x(n) = d_prime(n);
for i = n-1:-1:1
    x(i) = d_prime(i) - c_prime(i) * x(i+1);
end

end

function [L, U] = lu_banded(A, p, q)
% LU_BANDED 带状矩阵的LU分解
% 输入:
%   A - 带状矩阵
%   p - 下带宽
%   q - 上带宽

n = size(A, 1);
L = eye(n);
U = A;

for k = 1:n-1
    % 只处理带状范围内的元素
    i_start = max(1, k+1);
    i_end = min(n, k+p);
    
    for i = i_start:i_end
        if abs(U(k, k)) < eps
            error('主元为零');
        end
        L(i, k) = U(i, k) / U(k, k);
        
        j_start = max(1, k);
        j_end = min(n, k+q);
        U(i, j_start:j_end) = U(i, j_start:j_end) - L(i, k) * U(k, j_start:j_end);
    end
end

end

function x = solve_cyclic_tridiagonal(a, b, c, alpha, beta, d)
% SOLVE_CYCLIC_TRIDIAGONAL 求解循环三对角系统
% 使用Sherman-Morrison公式

n = length(b);

% 求解两个三对角系统
% 1) T*y = d
y = thomas_algorithm(a, b, c, d);

% 2) T*q = u，其中 u = [beta; 0; ...; 0; alpha]
u = zeros(n, 1);
u(1) = beta;
u(n) = alpha;
q = thomas_algorithm(a, b, c, u);

% Sherman-Morrison公式
% x = y - (v^T * y / (1 + v^T * q)) * q
% 其中 v = [1; 0; ...; 0; 1]
v_dot_y = y(1) + y(n);
v_dot_q = q(1) + q(n);

x = y - (v_dot_y / (1 + v_dot_q)) * q;

end